Development of VUMAT and VUHARD Subroutines for Simulating the Dynamic Mechanical Properties of Additively Manufactured Parts

Numerical modelling and simulation can be useful tools in qualification of additive manufactured parts for use in demanding structural applications. The use of these tools in predicting the mechanical properties and field performance of additive manufactured parts can be of great advantage. Modelling and simulation of non-linear material behaviour requires development and implementation of constitutive models in finite element analysis software. This paper documents the implementation and verification process of a microstructure-variable based model for DMLS Ti6Al4V (ELI) in two separate ABAQUS/Explicit subroutines, VUMAT and VUHARD, available for defining the yield surface and plastic deformation of materials. The verification process of the implemented subroutines was conducted for single and multiple element tests with varying prescribed loading conditions. The simulation results obtained were then compared with the analytical solutions at the same conditions of strain rates and temperatures. This comparison showed that both developed subroutines were accurate in predicting the flow stress of various forms of DMLS Ti6Al4V (ELI) under different conditions of strain rates and temperatures.


Introduction
Additive manufacturing (AM) technologies make it possible to design and fabricate lightweight metallic structural parts in real time. This is plausible, as components from 3D CAD models can be printed or produced directly using an electron beam source for electron beam melting (EBM) [1] or a laser source for selective laser melting (SLM)/direct metal laser sintering (DMLS) [1] on a powder bed table. Structural parts with simple to complex geometry in 3D CAD models can be fabricated via AM, since the technology offers great flexibility in design and manufacturing. Complex components with hollows and undercuts, such as turbine blades with internal cooling channels, are now being produced by AM [2]. Other complex components produced by this technology are patient-specific bio-implants such as dental prostheses and orthopaedic implants [3]. Therefore, the AM-produced components possess increased applications in the biomedical and aerospace industries.
As AM revolutionises manufacturing of functional components, much attention has focussed not only on accumulating exhaustive knowledge on the process itself, but also on resultant microstructures and mechanical properties of the manufactured parts. Steel, as the most common engineering structural material, is manifestly of high interest for AM. The grades of steel presently being produced by EBM and DMLS processes include maraging steel (18Ni300), stainless steels (17-4 PH, 15-5 PH, 321, 347, SAE 316L), and AISI 4340 steel [4]. The number of aluminium alloys produced by AM is still limited. The possible reason for this is that it is easy to machine aluminium, and furthermore the cost Numerical modelling and simulation can be useful to predict mechanical properties of materials, such as yield stresses and flow stresses, in those circumstances where actual experiments are too costly or even difficult to perform [18,19].
Various researchers have attempted to model and simulate the AM process and resulting microstructures. The finite element modelling (FEM) of an AM process to optimise the process parameters and therefore minimise residual stresses and distortion was reported in [20]. In this research work, the FEM software ABAQUS was employed to simulate an AM process, as it provides an interface allowing the user to define input data, such as element and heat input, that are a function of both position and time to achieve the process simulation of 3D parts. Marques et al. [21] used numerical analysis of the SLM process to assess the effects of laser scanning strategies on the residual stresses generated in the fabrication of Ti6Al4V parts. Chen et al. [22] used a 3D finite element model to simulate thermal behaviour in the SLM process and further investigated the effects of laser beam power and laser scanning speed on the thermal behaviour and solidification of SLM-produced high strength tool steel. Using the cellular automation (CA) method to describe grain growth, Zinoviev et al. [23] developed a two-dimensional numerical model to simulate the evolution of grain structure observed during the laser additive manufacturing process. Yang et al. [24] used a cellular automation model to simulate microstructural solidification and solid-state phase transformation of the Ti6Al4V alloy under various spatially variable thermal cycles of the SLM process. In their work, the morphology and size of β-grains and martensite simulated by the model agreed well with experimental results. Artificial neutral network (ANN) models have also been utilised to optimise AM process parameters and predict the mechanical behaviour of 3D-printed parts [25,26]. These types of models depend on huge amounts of experimental data for model-training purposes, and constitutive mathematical expressions describing the material properties are not required. This last bit makes it difficult to use ANN models to develop numerical models for predicting the mechanical properties of AM parts.
This study focuses on developing a numerical model based on the dominant mechanisms of deformation of microstructures of AM Ti6Al4V (ELI), over a broad-spectrum of temperatures and strain rates. To carry out modelling and simulation of a non-linear analysis considering large displacements and elastoplastic material behaviour, user-defined material models are required. This is so in the present case where in-built models in finite element analysis (FEA) packages are inadequate. The commercial package of ABAQUS provides user interfaces via a FORTRAN compiler linking to the main program to allow the user to formulate and incorporate user-defined material models [27].
Therefore, the developed, calibrated, and validated microstructure-and dislocationbased analytical constitutive model reported in the authors' previous work in [28,29] is implemented in this study as a user material subroutine in ABAQUS using VUMAT and VUHARD subroutines. The VUHARD subroutine is used to describe the plastic deformation of a material, while the VUMAT subroutine describes the complete (elasticplastic) behaviour of a material. Both subroutines use an explicit integration scheme (forward Euler method), which calculates the state of the variables from the model at the current time [27]. Implementation of the VUHARD subroutine is easy and straightforward and requires only the constitutive model and its derivatives. Implementation of the VUMAT subroutine, on the other hand, is involved and requires a clear understanding of the von Mises yield criterion and the radial return algorithm method used to evaluate the increment of state variables. A summary of the von Mises yield criterion and the radial return algorithm used in most FEA programs is presented in Appendix A.
The codes developed in the present work for each of the two ABAQUS user subroutines are presented here as Supplementary Material. The implemented subroutines were verified using single and multiple element tests with varying prescribed loading. The simulation results obtained were then compared with the analytical solutions developed and reported elsewhere by the author [28] and showed good accuracy for both the developed VUMAT and VUHARD subroutines applied to high-strain-rate numerical simulations.

A Microstructure and Dislocation-Based Model for DMLS Ti6Al4V (ELI)
This section presents an analytical model that uses the critical microstructural parameters of various forms of DMLS Ti6Al4V (ELI) to predict the mechanical properties of the alloy over a wide range of strain rate and temperature. Details on formulation of the model were presented in the authors' previous work published in [29]. The total flow stress of a given microstructure of DMLS Ti6Al4V (ELI) is dependent on internal and external state variables. The internal state variables include the average grain size and the density of dislocations, while the external state variables are field variables of strain, strain rate, and temperature. The incorporation of mechanical threshold stress (MTS), the Hall-Petch relationship, the Taylor strain hardening law, and viscous drag effects at a high strain rate were used in [29] to describe the yielding and flow stress of DMLS Ti6Al4V (ELI) microstructures and a wide range of field variables. The total flow stress developed this way is of the following form: The first part of the foregoing expression accounts for the MTS where the symbol σ o is the mechanical threshold stress or the value of the thermal stress at 0 K, µ is the shear modulus dependent on temperature T, and µ o is the shear modulus at absolute zero temperature (0 K). The symbols k b , b, g oi , .
ε o , and . ε are the Boltzmann's constant, the Burgers vector, the material constant, the reference strain rate, and the testing strain rate, respectively. The constants p and q are model fitting parameters. The second part of Equation (1) describes the upsurge of flow stress as the plastic strain rate ( . ε p ) increases because of viscous drag effects, and the symbols ζ and χ in it are viscous drag stress calibration parameters. The third part of this equation represents the inverse relationship between the average thickness of α-laths, t, and the yield stress of a material, where K H−P is the Hall-Petch constant. The fourth and last part of this equation is the integrated Taylor hardening law which represents the accumulation of dislocations with increase in plastic strain. In this part of the equation, the symbol α is a dimensionless parameter of magnitude ranging from 0.2 to 0.4 for different materials. The parameter M is the Taylor factor which relates the shear flow stress τ of a single crystal to the uniaxial flow stress σ of a polycrystal. The constants h and k 2 stand for the coefficient of dislocation accumulation and annihilation, respectively, while ρ 0 is the total initial dislocation density in a particular material's microstructure.

Materials
The calibration and validation of the flow stress model in Equation (1) were reported in [29] and required material properties obtained from analysis of microstructures and experimental tests. Cylindrical rods with a height of 80 mm and diameter 6 mm were produced from Ti6Al4V (ELI) (grade 23) powder by the DMLS process using an EOS M280 machine. The chemical composition (in wt.%) of this powder as provided by the suppliers, TLS Technik GmbH (Bitterfeld-Wolfen, Germany), is Al: 6.34, V: 3.944, Fe: 0.25, O: 0.082, C: 0.006, N: 0.006, H: 0.00, and the balance Ti.
These rods were first stress relieved at a temperature of 650 • C for 3 h and then furnace cooled to room temperature. The stress relieved rods were then subdivided into three groups designated hereinafter as samples C, D, and E for different heat treatment processes. The samples designated C were heat-treated at 800 • C for 2.5 h and then furnace-cooled to room temperature. The samples designated D were duplex-annealed at 940 • C for 2 h and then furnace-cooled, followed by heat treatment at 750 • C for 2 h and then finally furnace-cooled. The samples designated E were heat treated above the β-transformation temperature at 1020 • C for 2.5 h and then furnace-cooled to room temperature. The optical These rods were first stress relieved at a temperature of 650 °C for 3 h and then furnace cooled to room temperature. The stress relieved rods were then subdivided into three groups designated hereinafter as samples C, D, and E for different heat treatment processes. The samples designated C were heat-treated at 800 °C for 2.5 h and then furnacecooled to room temperature. The samples designated D were duplex-annealed at 940 °C for 2 h and then furnace-cooled, followed by heat treatment at 750 °C for 2 h and then finally furnace-cooled. The samples designated E were heat treated above the β-transformation temperature at 1020 °C for 2.5 h and then furnace-cooled to room temperature. The optical images and secondary electron images of the microstructures obtained after these different heat treatment strategies are shown in Figure 1a-f. The detailed description of these microstructures can be found in [10,28,30]. The average α-grain size in each of these microstructures was determined by the line intercept method and results presented in [10]. The total initial dislocation density in these heattreated samples was also determined by X-ray diffraction peak broadening analysis and The detailed description of these microstructures can be found in [10,28,30]. The average α-grain size in each of these microstructures was determined by the line intercept method and results presented in [10]. The total initial dislocation density in these heattreated samples was also determined by X-ray diffraction peak broadening analysis and the results reported in [30]. As seen in Equation (1), the average grain size and the initial dislocation density are the two critical internal state variables in any given microstructure of DMLS Ti6Al4V (ELI) required to describe the flow stress. These important numerical model input parameters for different forms of DMLS Ti6Al4V (ELI) are shown in Table 1.

Experimental Tests
The flow stress curves of samples C, D, and E over a wide range of strain rates and temperatures were obtained using the Kolsky bar commonly known as the split-Hopkinson pressure bar (SHPB). The tests were conducted at strain rates of 750, 1500, and 2450 s −1 and the tests at each strain rate carried out at the three different temperatures of 25, 200, and 500 • C. The results obtained in these tests were reported in [28], while the calibration and validation of the microstructural sensitive model presented in Section 2 using flow curves in [28] were reported in [29]. A summary of the model-calibrated fitting parameters obtained from this analytical modelling and other physical constants used and obtained from literature is presented in Table 2. Implementation of the VUHARD user subroutine for the microstructure-based constitutive formulation of flow stress was done by first obtaining the derivatives of this constitutive model with respect to appropriate variables. In this constitutive model, presented here as Equation (1), the total flow stress is seen to be dependent on three external state variables (ε p , . ε p , and T), whilst the internal state variables can be obtained from the microstructure. The temperature-dependence of the constitutive model was taken care of by the decrease in shear modulus (µ) with an increase in temperature using the equation for shear modulus shown in Table 2.
Therefore, in this study, implementation of the VUHARD subroutine for the constitutive model required only two analytical derivatives of D 1 = ∂σ/∂ ε p and D 2 = ∂σ/∂ . ε p , which were then evaluated numerically to allow quadratic convergence. These derivatives represented the change in the calculated total flow stress for the inputs of equivalent plastic strain rate and equivalent plastic strain. The term "equivalent" is used here to denote the scalar quantity of strain and strain rate obtained from the transformation of the strain tensor using the von Mises criterion (presented and discussed in Appendix A), which is the Materials 2022, 15, 372 7 of 32 common criterion used in ABAQUS. The partial derivative ∂σ n+1 /∂ε n+1 p was determined numerically from the following equation: Similarly, the other derivative, with respect to plastic strain rate, . ε n+1 p , was evaluated from the following equation: The VUHARD user subroutine written in this work is included in this paper as Supplementary Material in S1. The coding for this subroutine was generated using FORTRAN 90 and saved in a .for file, formatted to be compatible with ABAQUS 2020. The first algorithm for the subroutine was used to generate the user material properties described as user material constants. A total of 17 user material constants, derived from Equation (1) and listed in Table 3, were generated. These excluded the material density, elastic modulus, and Poisson's ratio that were provided separately in ABAQUS/CAE. The flow stress component that is dependent on the average thickness of α-laths (t) was computed based on the Hall-Petch relationship and was input into the subroutine as an average variable of α-lath thickness (a variable that depends on the form of the alloy).
The previous and current values of equivalent plastic strain in a step were defined as old and new solution-dependent variables (SDV1), respectively. Thus, the increment of plastic strain was obtained as the difference between these two states. The equivalent plastic strain rate component of Equations (1) and (3) was expressed from the increment in equivalent plastic strain (∆ε p ) for a time increment size (∆t) as follows: Equations (1)-(3) were evaluated at given material points and arrays containing values of equivalent yield stress, plastic strain, and strain rate, generated at related material points.

Implementation of the Microstructure-Based Constitutive Flow Stress Model as a VUMAT Subroutine
To make it easier to implement the microstructure-and dislocation-based constitutive model as a VUMAT user material subroutine, the implementation algorithm was written following the logic of the flow chart shown in Figure 2. Coding for this subroutine was based on the theory of plasticity and radial return mapping presented in Appendix A.  If , then the state response was elastic, and the radial return algorithm that is inherent in the subroutine to correct the plastic response was skipped. The stress tensor and related deviatoric stress remain unchanged at any value of q. Therefore, the last step of the algorithm to compute internal energy and dissipation of inelastic en- Similar to VUHARD, the coding for VUMAT was generated in FORTRAN 90 and saved as a .for file for execution in ABAQUS 2020. The codes written are presented in S2 in the Supplementary File. The second block in the algorithm was used to obtain the model constants of the material. Both elastic and plastic properties are required in this subroutine's main program, which is different from VUHARD where only the plastic properties are provided in the main program and the elastic properties are provided in ABAQUS/CAE. Thus, the 17 material constants from Equation (1) in addition to the material's elastic modulus (E) and Poisson's ratio (ν) were used for this subroutine. These material constants are presented in Table 3. The two Lamé's constants (µ and λ) were computed from the elastic modulus (E) and Poisson's ratio (ν) using the following two expressions for isotropic hardening: The principal components of the stress tensor were first obtained at the reference time t + ∆t = 0 from the following trace of principal strain increments: At time (t + ∆t > 0), the VUMAT algorithm was split in three parts. In the first part, the values of equivalent plastic strain (ε n p ), equivalent plastic strain rate ( . ε n p ), and yield stress σ n y (given by the constitutive model at ε n p and . ε n p and at a prescribed test temperature (T)) at the beginning of the increments at each integration point were defined as SDV1, SDV2, and SDV3, respectively. The values of shear modulus used in the constitutive model to define the flow stress were computed at prescribed test temperatures before obtaining the value of flow stress in this part.
In the second part, the theory of plasticity presented in Appendix A was used to compute equivalent von Mises trial stresses q. The third part aimed at comparing the values of equivalent von Mises trial stresses obtained with the values of yield stress (from the constitutive model given by Equation (1)) obtained at the beginning of each increment such that: If q < σ n y , then the state response was elastic, and the radial return algorithm that is inherent in the subroutine to correct the plastic response was skipped. The stress tensor and related deviatoric stress remain unchanged at any value of q. Therefore, the last step of the algorithm to compute internal energy and dissipation of inelastic energy was initiated. It is important to note that internal energy, abbreviated as ALLIE in ABAQUS, is the sum of the recoverable elastic strain energy and the energy dissipated through plastic deformation [27]. Since there was no plastic deformation, the inelastic energy remained zero, and the internal energy only included the elastic strain energy component. If q ≥ σ n y , then the stress state was outside the yield surface and was plastic. Thus, the radial return mapping previously described in Appendix A was initiated to compute the equivalent plastic strain and therefore return the predicted stress back to the expanded yield surface of the material. At this stage, strain hardening was expected since the loading was beyond the elastic limit.
The final part of this algorithm updated the SDVs and computed the new internal energy and inelastic dissipated energy. The updated SDV1, SDV2, and SDV3 at each stage of computation, which were now the new values of equivalent plastic strain ε n+1 p , the new equivalent plastic strain rate . ε n+1 p , and the new yield stress σ n+1 y , were stored and then used in the next increment as current values. For linear isotropic materials undergoing small strains, the strain energy is of the following form [27]: From Equation (7), the total new specific internal energy (inter.E new ) and specific inelastic dissipation energy (inel.E new ) were expressed as

Results and Discussions
The microstructure-and dislocation-based constitutive model implemented as VUHARD and VUMAT subroutines was verified by determining the closeness of the results of simulation using ABAQUS to the analytical solutions based on Equation (1). The user material constants in this equation were those that were used to fit the equation to experimental results for samples C, D, and E, which have been presented in Tables 1 and 2. Various benchmark tests were used here ranging from single element tests to multiple element tests.

The Single Element Tests
The verification process adopted in the present research was to first test the material model developed in ABAQUS using single eight-node linear brick element tests. This is an easy and practical technique to examine the sensitivity and accuracy of an element to external loading. Here, the single eight-node linear brick element tests were performed using the prescribed velocity loading conditions both in compression and tension, as shown in Figure 3a,b, respectively.  The roller support condition to restrain the surface in a direction parallel to the loading direction was applied to the four nodes of the surface opposite to the surface onto which load was applied. The dimensions of the cube were arbitrarily taken as (10 × 10 × 10) mm and the cube meshed using the eight-node linear brick C3D8R element shown in Figure 3c, with reduced integration and hourglass control settings. This type of mesh is simple and appropriate for the cubic element models of components with regular straightedged shapes.
The reduced integration elements normally have a single integration point located at the element's centroid rather than full integrated element types and hence give rise to The roller support condition to restrain the surface in a direction parallel to the loading direction was applied to the four nodes of the surface opposite to the surface onto which load was applied. The dimensions of the cube were arbitrarily taken as (10 × 10 × 10) mm and the cube meshed using the eight-node linear brick C3D8R element shown in Figure 3c, with reduced integration and hourglass control settings. This type of mesh is simple and appropriate for the cubic element models of components with regular straightedged shapes.
The reduced integration elements normally have a single integration point located at the element's centroid rather than full integrated element types and hence give rise to reduced computation time. However, reduced integration elements tend to be too flexible, since they suffer from their own numerical problem called hourglassing. Hourglassing arises where elements are unable to resist bending and hence cannot store strain energy for bending deformation. In ABAQUS hourglass control, a small amount of artificial stiffness is introduced in elements to limit the propagation of the hourglass phenomenon [27].
Three      With dynamic explicit integration of instantaneous imposed velocities, step times are important since long step times could cause excessive deformation of elements. In some instances, the opposite faces may cross over each other (in compression) causing the job to exit with an error, whilst for very small step times, the elements may not reach the plastic deformation state. For single element tests, the dynamic explicit integration option with total step times (t) of 0.03 s, 0.003 s, and 0.001 s for velocities of 0.1 m/s, 1 m/s, and 4 m/s, respectively, were found to give sufficient deformation without excessive distortion of the element. For velocities beyond 4 m/s for this kind of numerical model, the simulation aborted even for very small step times due to excessive deformation of the model. Figure 6 shows the change in the values of equivalent von Mises stresses, plastic strains, and plastic strain rates for samples C with increasing simulation time in VUHARD and VUMAT for a velocity of 4 m/s. In both subroutines, the plastic strain is initially zero (during elastic deformation) before increasing uniformly to a maximum value, and therefore the rate of deformation (plastic strain rate) is constant during this period, as seen in Figure 6. Both VUMAT and VUHARD show fairly similar values of average plastic strain rate of approximately 350 s −1 in the figure. The equivalent plastic strain and flow stress curves generated from the two subroutines also coincide for the entire and a large part of the simulation step time, respectively. of the element. For velocities beyond 4 m/s for this kind of numerical model, the simulation aborted even for very small step times due to excessive deformation of the model. Figure 6 shows the change in the values of equivalent von Mises stresses, plastic strains, and plastic strain rates for samples C with increasing simulation time in VUHARD and VUMAT for a velocity of 4 m/s. In both subroutines, the plastic strain is initially zero (during elastic deformation) before increasing uniformly to a maximum value, and therefore the rate of deformation (plastic strain rate) is constant during this period, as seen in Figure  6. Both VUMAT and VUHARD show fairly similar values of average plastic strain rate of approximately 350 s −1 in the figure. The equivalent plastic strain and flow stress curves generated from the two subroutines also coincide for the entire and a large part of the simulation step time, respectively. As discussed in the authors' previous work reported in [29], there are four critical parameters of the constitutive model implemented in the subroutines in this research that influence the shape of the flow stress curve, namely ℎ, , , and . The sensitivity of these critical parameters in the numerical models for samples C developed here was investigated at a velocity of 4 m/s, and the results were compared with those of the analytical solutions at the same conditions for both cases. These results are presented in Figures 7 and 8. As discussed in the authors' previous work reported in [29], there are four critical parameters of the constitutive model implemented in the subroutines in this research that influence the shape of the flow stress curve, namely h, k 2 , ρ o , and T. The sensitivity of these critical parameters in the numerical models for samples C developed here was investigated at a velocity of 4 m/s, and the results were compared with those of the analytical solutions at the same conditions for both cases. These results are presented in Figures 7 and 8.
As seen in these figures, the results from the numerical models derived in VUMAT and VUHARD, and those from the analytical solutions developed here, coincide for a large part of the flow curves. The opposing effect of parameters h and k 2 is well articulated in the numerical model where the flow stress increases with decreasing values of the parameter k 2 , whereas it increases with increasing values of the parameter h, as seen in Figure 7. Similar observations were reported in [29]. Figure 8 shows the effects of initial dislocation density and temperature on the numerical solutions obtained from the VUMAT and VUHARD subroutines developed here and the analytical solutions, all at the same strain rate of 350 s −1 . It is clear from Figure 8a that for very high initial dislocation densities, the flow stress curve is characterised by high initial peak stresses followed by a decrease in flow stress. For very low initial dislocation densities, the peak stresses are lower, and the material experiences strain hardening. Similar observations were noted in [29], and the theory behind this phenomenon was presented.
The effects of simulation temperatures on the flow stress curve are very conspicuous in Figure 8b, where the values of flow stress are seen to decrease with increase in temperature. This is also consistent with observations made in [29]. The preceding discussion and the observations made from Figures 8 and 9 show good sensitivity of the numerical models developed here in VUMAT and VUHARD subroutines to the critical parameters that dictate the shape of the flow stress curve in the DMLS Ti6Al4V (ELI) alloy. The close coincidence of the curves for the numerical models with those for the analytical models is also noted. This creates confidence in the use of these models in the simulation of flow stresses for sample types C, D, and E at various imposed velocities.  As seen in these figures, the results from the numerical models derived in VUMAT and VUHARD, and those from the analytical solutions developed here, coincide for a large part of the flow curves. The opposing effect of parameters ℎ and is well articulated in the numerical model where the flow stress increases with decreasing values of the parameter , whereas it increases with increasing values of the parameter ℎ, as seen in Figure 7. Similar observations were reported in [29]. Figure 8 shows the effects of initial dislocation density and temperature on the numerical solutions obtained from the VUMAT and VUHARD subroutines developed here and the analytical solutions, all at the same strain rate of 350 s −1 . It is clear from Figure 8a that for very high initial dislocation densities, the flow stress curve is characterised by high initial peak stresses followed by a decrease in flow stress. For very low initial dislocation densities, the peak stresses are As seen in these three figures, the scaling down and upturn of flow stress at high temperature and with an increase in strain rate, respectively, is well articulated. For the two cases of strain rate and temperature, the results from the numerical models based on the VUMAT and VUHARD subroutines and those derived from the analytical solution developed here are nearly overlapping over most of the equivalent plastic strain range. This confirms that the two subroutines work well in the loading conditions applied.
As seen in these three figures, the scaling down and upturn of flow stress at high temperature and with an increase in strain rate, respectively, is well articulated. For the two cases of strain rate and temperature, the results from the numerical models based on the VUMAT and VUHARD subroutines and those derived from the analytical solution developed here are nearly overlapping over most of the equivalent plastic strain range. This confirms that the two subroutines work well in the loading conditions applied.  Both the compression and tensile test simulations showed the same values of flow stress at all strain rates and temperatures, which is anticipated considering the criterion for isotropic yielding used to develop and implement the constitutive models. Even though this was the case, for the same simulation time, there was a significant difference in the maximum equivalent plastic strains attained in the two load conditions, with values of about 0.35 in compression and about 0.25 in tension. The elastic moduli in tension and compression for metals and alloys are normally taken to be the same. Additionally, failure due to the formation of adiabatic shear bands (ASBs) occurs in both high-strain-rate compression and tension. Failure due to monotonic tensile stress occurs first through the formation of pores and cracks and their horizontal spread in the initial brittle fracture and, thereafter, necking before final tearing along the ASBs. Compressive stresses, on the other hand, cause materials to reduce in length and increase in cross section, and failure in most cases is because of the initiation of ASBs. These distinct deformation mechanisms in compression and tension suggest different values of plastic strain for the same applied load, as seen in Figures 9 and 11.  The necking failure mechanism in tension can be studied better using multiple elements, since multiple element models consist of several integration points, and therefore results from various sections of the model can be obtained. A specific geometry of the model can also be used to allow early necking.

Multiple Elements Tests
The necking of a circular bar test is an example commonly used in finite element analysis to investigate the necking and softening of a tensile bar [31]. The test was used in the present research to investigate the performance of the VUHARD and VUMAT subroutines developed here under the conditions of necking in tensile loading. The numerical model used was that of the asymmetrical cylinder shown in Figure 12a, with large and small diameters (on the fixed end and at the surface of application of load) of 10 mm and 8 mm, respectively, and a height of 10 mm. The roller boundary condition with y-axis rotational symmetry was applied on the side of the model with a large diameter to constrain the displacement only in the y-direction and to allow rotation only about the y-axis. Tensile loading of the cylinder was realised through imposed instantaneous velocities of 0.1 m/s, Both the compression and tensile test simulations showed the same values of flow stress at all strain rates and temperatures, which is anticipated considering the criterion for isotropic yielding used to develop and implement the constitutive models. Even though this was the case, for the same simulation time, there was a significant difference in the maximum equivalent plastic strains attained in the two load conditions, with values of about 0.35 in compression and about 0.25 in tension. The elastic moduli in tension and compression for metals and alloys are normally taken to be the same. Additionally, failure due to the formation of adiabatic shear bands (ASBs) occurs in both high-strain-rate compression and tension. Failure due to monotonic tensile stress occurs first through the formation of pores and cracks and their horizontal spread in the initial brittle fracture and, thereafter, necking before final tearing along the ASBs. Compressive stresses, on the other hand, cause materials to reduce in length and increase in cross section, and failure in most cases is because of the initiation of ASBs. These distinct deformation mechanisms in compression and tension suggest different values of plastic strain for the same applied load, as seen in Figures 9 and 11.
The necking failure mechanism in tension can be studied better using multiple elements, since multiple element models consist of several integration points, and therefore results from various sections of the model can be obtained. A specific geometry of the model can also be used to allow early necking.

Multiple Elements Tests
The necking of a circular bar test is an example commonly used in finite element analysis to investigate the necking and softening of a tensile bar [31]. The test was used in the present research to investigate the performance of the VUHARD and VUMAT subroutines developed here under the conditions of necking in tensile loading. The numerical model used was that of the asymmetrical cylinder shown in Figure 12a, with large and small diameters (on the fixed end and at the surface of application of load) of 10 mm and 8 mm, respectively, and a height of 10 mm. The roller boundary condition with -axis rotational symmetry was applied on the side of the model with a large diameter to constrain the displacement only in the -direction and to allow rotation only about theaxis. Tensile loading of the cylinder was realised through imposed instantaneous velocities of 0.1 m/s, 1 m/s, and 4 m/s along the -axis (the axial direction of the cylinder), shown in Figure 12b, each at a time. These four different sizes of meshes yielded 475, 3550, 28,320, and 126,038 linear hexahedral elements of C3D8R type, in order from the largest mesh size to the smallest. Figures 13-16 show the distribution of equivalent von Mises stresses and plastic strains obtained using the VUHARD and VUMAT subroutines developed here, for the four different mesh sizes. These four different sizes of meshes yielded 475, 3550, 28,320, and 126,038 linear hexahedral elements of C3D8R type, in order from the largest mesh size to the smallest. Figures 13-16 show the distribution of equivalent von Mises stresses and plastic strains obtained using the VUHARD and VUMAT subroutines developed here, for the four different mesh sizes.          A summary of computation times for simulation of the necking bar in this set of simulations is shown in Table 4. The VUMAT subroutine is seen in this table to have more computation time and increments in comparison to VUHARD. The difference in computation time and increments between these two subroutines is possibly due to that fact that they use different integration schemes (discussed later in this section) to update the state variables that could influence the stable time increment. The computation time is also seen to increase with mesh density. For all four sizes of mesh, the VUHARD and VUMAT subroutines display similar distributions of equivalent plastic strain (ε p ) and von Mises stress (σ p ). From the legends in (a) and (c) of Figure 13 to Figure 16, all four numerical models are seen to experience plastic deformation. The maximum values of equivalent plastic strain and equivalent von Mises stress are located at the tip labelled y 2 and decline gradually towards y 1 where they are minimum, as seen in the contours in these figures.
It should be noted that the shape of the undeformed model shown in Figure 12 is that of a truncated cone with the cross-section areas that decrease linearly towards the top surface. Thus, for a given load that is applied at the top surface of the model, the values of equivalent von Mises stress and strain will decline as the cross section increases towards the bottom surface of the original cone.
To better understand the distribution of these two variables for different mesh densities, their values were obtained from elements along the surface from y 1 to y 2 at the end of each simulation. These values were then plotted against the normalised length of the deformed model. The resulting curves for different mesh densities are shown in Figure 17, while the results of mesh convergence studies on various sections of the model are presented in Figure 18.  Table 4. The VUMAT subroutine is seen in this table to have more computation time and increments in comparison to VUHARD. The difference in computation time and increments between these two subroutines is possibly due to that fact that they use different integration schemes (discussed later in this section) to update the state variables that could influence the stable time increment. The computation time is also seen to increase with mesh density. For all four sizes of mesh, the VUHARD and VUMAT subroutines display similar distributions of equivalent plastic strain ( ) and von Mises stress ( ). where they are minimum, as seen in the contours in these figures. It should be noted that the shape of the undeformed model shown in Figure 12 is that of a truncated cone with the cross-section areas that decrease linearly towards the top surface. Thus, for a given load that is applied at the top surface of the model, the values of equivalent von Mises stress and strain will decline as the cross section increases towards the bottom surface of the original cone.
To better understand the distribution of these two variables for different mesh densities, their values were obtained from elements along the surface from to at the end of each simulation. These values were then plotted against the normalised length of the deformed model. The resulting curves for different mesh densities are shown in Figure  17, while the results of mesh convergence studies on various sections of the model are presented in Figure 18.   Both the equivalent plastic strains ( ) and von Mises stresses ( ) show distinct profiles for the three different mesh densities of the model in Figure 17a,b. The profiles exhibit three typical zones, as shown in Figure 17a, which differ in range for different mesh sizes. The first zone shows a slow increase in at increasing rates and its ranges from normalised lengths of 0-0.26, 0-0.47, 0-0.60, and 0-0.60 for models with 2 mm, 1 mm, 0.5 mm, and 0.25 mm mesh sizes, respectively. The slow increase of equivalent strain in this zone is attributed to the large cross-section areas of the model and, therefore, low equivalent von Mises stress acting on elements in the region. The second zone shows a higher increase in with points of inflexion. It ranges from normalised lengths of 0.26 -0.63, 0.47-0.70, 0.60-0.73, and 0.60-0.73 in order of decreasing mesh size. This is the zone that precedes necking. A rapid decrease in the cross-sectional area around this zone resulted in a rapid increase in equivalent von Mises stress and strain. This stage forms the shoulder of the "neck". The third zone is the necking section of the model and shows an increase in at constant rates, and it ranges from normalised lengths of 0.47-1, 0.70-1, 0.73-1, and 0.73-1 in order of increasing mesh density.
The mesh convergence analysis shown in Figure 18a for the values of equivalent plastic strain obtained at various sections of the model gives rise to several observations. There is convergence of the values of equivalent plastic strain as the mesh density increases in the second and third zones. Lower and higher values of equivalent plastic strain from the VUHARD subroutine in the second and third zones, respectively, are also noted. There is insignificant variation of equivalent plastic strain as the mesh density increases in the first zone, and the results from the VUHARD and VUMAT subroutines are perfectly superimposed in this zone.
From Figure 17b, the influence of mesh density on the distribution of equivalent von Mises stresses is seen to be significant between the normalised lengths of 0 and 0.14. The values of equivalent von Mises stresses are seen to increase with increasing mesh size in this range for both the results obtained using the VUMAT and VUHARD subroutines. However, with 0.5 mm and 0.25 mm mesh sizes, the results from the two subroutines overlap. In the zone between normalised lengths of 0.14 and 0.57, there is a small spread of results obtained for different mesh densities. The variation between the values of equivalent von Mises stresses obtained with coarse and refined mesh models in this region is less than 0.1% of the former. Similar to the preceding zone, the 0.5 mm and 0.25 mm mesh sizes in this zone give results that are indistinguishable. The results obtained in the final zone for normalised lengths of 0.57 to 1 are seen in Figure 17b to be superimposed on top Both the equivalent plastic strains (ε p ) and von Mises stresses (σ p ) show distinct profiles for the three different mesh densities of the model in Figure 17a,b. The ε p profiles exhibit three typical zones, as shown in Figure 17a, which differ in range for different mesh sizes. The first zone shows a slow increase in ε p at increasing rates and its ranges from normalised lengths of 0-0.26, 0-0.47, 0-0.60, and 0-0.60 for models with 2 mm, 1 mm, 0.5 mm, and 0.25 mm mesh sizes, respectively. The slow increase of equivalent strain in this zone is attributed to the large cross-section areas of the model and, therefore, low equivalent von Mises stress acting on elements in the region. The second zone shows a higher increase in ε p with points of inflexion. It ranges from normalised lengths of 0.26 -0.63, 0.47-0.70, 0.60-0.73, and 0.60-0.73 in order of decreasing mesh size. This is the zone that precedes necking. A rapid decrease in the cross-sectional area around this zone resulted in a rapid increase in equivalent von Mises stress and strain. This stage forms the shoulder of the "neck". The third zone is the necking section of the model and shows an increase in ε p at constant rates, and it ranges from normalised lengths of 0.47-1, 0.70-1, 0.73-1, and 0.73-1 in order of increasing mesh density.
The mesh convergence analysis shown in Figure 18a for the values of equivalent plastic strain obtained at various sections of the model gives rise to several observations. There is convergence of the values of equivalent plastic strain as the mesh density increases in the second and third zones. Lower and higher values of equivalent plastic strain from the VUHARD subroutine in the second and third zones, respectively, are also noted. There is insignificant variation of equivalent plastic strain as the mesh density increases in the first zone, and the results from the VUHARD and VUMAT subroutines are perfectly superimposed in this zone.
From Figure 17b, the influence of mesh density on the distribution of equivalent von Mises stresses is seen to be significant between the normalised lengths of 0 and 0.14. The values of equivalent von Mises stresses are seen to increase with increasing mesh size in this range for both the results obtained using the VUMAT and VUHARD subroutines. However, with 0.5 mm and 0.25 mm mesh sizes, the results from the two subroutines overlap. In the zone between normalised lengths of 0.14 and 0.57, there is a small spread of results obtained for different mesh densities. The variation between the values of equivalent von Mises stresses obtained with coarse and refined mesh models in this region is less than 0.1% of the former. Similar to the preceding zone, the 0.5 mm and 0.25 mm mesh sizes in this zone give results that are indistinguishable. The results obtained in the final zone for normalised lengths of 0.57 to 1 are seen in Figure 17b to be superimposed on top of one another. It is important to note here that for large values of equivalent plastic strain, the flow stress given by the constitutive model developed in this study saturates (no strain hardening), which elucidates the convergence of all these curves in the final zone. This is clearly shown in Figure 18b at a normalised length of the model of 0.7.
The values of equivalent von Mises stresses obtained from the VUHARD subroutine are higher and lower than those from the VUMAT subroutine in the first and second zones at lower mesh densities, respectively. The reverse is true for higher mesh densities, with cross-over points at mesh densities of about 0.85 and 1.2 for the two zones, respectively, while in the last zone, the results from the two subroutines are indistinguishable. Although the results from these two subroutines are very close to one another, the small differences seen in the preceding two figures could be a result of the different integration schemes employed by these two subroutines. The VUHARD subroutine uses an explicit central difference time integration rule built into ABAQUS to generate plastic strain at Gauss points [32], while the VUMAT subroutine uses the radial return method, which belongs to backward or forward integration algorithm methods, in which equivalent plastic strains are evaluated at Gauss points. The difference in their time integration schemes suggests variations in the results obtained due to inherent integration errors associated with each integration scheme. It is important to note that the Euler and central difference methods of solving differential equations are more accurate if the step time is small. The error in the forward and backward integration algorithms is proportional to the size of the time step, while for the central difference method, the error is equal to the square of the size of the time step. Thus, for small time steps (<1 s), the central difference time integration scheme is more accurate than the forward and backward integration methods.
The plots of equivalent plastic strains against equivalent von Mises stresses obtained from the two subroutines are shown in Figure 19 for four different mesh sizes. As seen in this figure, the curves are nearly overlapping for the better part of the flow stress curves. The significant effects of mesh size are revealed at the initial values of equivalent plastic strain (circled in Figure 19), with the 2 mm mesh size giving a much higher value of about 0.083 equivalent plastic strain in comparison to the 0.25 mesh size, which has a value of about 0.02. However, the curves for 0.25 mm and 0.5 mm mesh sizes are indistinguishable at all values of strain, which is consistent with the observations made from Figures 17 and 18. of one another. It is important to note here that for large values of equivalent plastic strain, the flow stress given by the constitutive model developed in this study saturates (no strain hardening), which elucidates the convergence of all these curves in the final zone. This is clearly shown in Figure 18b at a normalised length of the model of 0.7. The values of equivalent von Mises stresses obtained from the VUHARD subroutine are higher and lower than those from the VUMAT subroutine in the first and second zones at lower mesh densities, respectively. The reverse is true for higher mesh densities, with cross-over points at mesh densities of about 0.85 and 1.2 for the two zones, respectively, while in the last zone, the results from the two subroutines are indistinguishable. Although the results from these two subroutines are very close to one another, the small differences seen in the preceding two figures could be a result of the different integration schemes employed by these two subroutines. The VUHARD subroutine uses an explicit central difference time integration rule built into ABAQUS to generate plastic strain at Gauss points [32], while the VUMAT subroutine uses the radial return method, which belongs to backward or forward integration algorithm methods, in which equivalent plastic strains are evaluated at Gauss points. The difference in their time integration schemes suggests variations in the results obtained due to inherent integration errors associated with each integration scheme. It is important to note that the Euler and central difference methods of solving differential equations are more accurate if the step time is small. The error in the forward and backward integration algorithms is proportional to the size of the time step, while for the central difference method, the error is equal to the square of the size of the time step. Thus, for small time steps (<1 s), the central difference time integration scheme is more accurate than the forward and backward integration methods.
The plots of equivalent plastic strains against equivalent von Mises stresses obtained from the two subroutines are shown in Figure 19 for four different mesh sizes. As seen in this figure, the curves are nearly overlapping for the better part of the flow stress curves. The significant effects of mesh size are revealed at the initial values of equivalent plastic strain (circled in Figure 19  It is not surprising to see the flow stress saturating at the equivalent plastic strain of 0.4 and onward for the VUMAT and VUHARD subroutines, since the implemented model was not developed to take care of unloading conditions or even evolution of damage. During the initial stage of plastic deformation, materials experience work-hardening where flow stress increases with increasing strain. Work-hardening normally occurs due to generation and accumulation of dislocations with increasing plastic strain. However, as previously reported in the authors' previous work in [28], materials loaded at a high strain rate experience "thermal softening". Thermal softening causes dynamic recovery where dislocations begin to annihilate. If the rate of annihilation of dislocations equals the rate of generation of dislocations by strain hardening, then the flow stress saturates. However, if the rate of annihilation of dislocations exceeds the rate at which dislocations are generated, materials experience decreasing flow stress with increasing plastic strain, a phenomenon commonly referred to as "strain softening". The threshold value of the accumulation of dislocations for the three forms of DMLS Ti6Al4V (ELI) samples was calibrated as 8.3 × 10 15 m −2 in [29]. The formation of ASBs is the precursor to material failure due to the formation and growth of microvoids and cracks along them [28]. The failure processes at high strain rates that lead to instabilities can be reflected in the numerical models through extension of the analytical model to predict the evolution of such damage. However, this was not within the scope of the present research, and future research is expected to focus on integrating the evolution of damage into the constitutive numerical model developed here.
The foregoing graphs and related discussions suggest that mesh density does influence the strain and stress distribution in a numerical model. However, this does not appear to significantly affect the VUMAT and VUHARD subroutines differently. This is evident from the fact that for the four different mesh densities used here, the VUMAT and VUHARD subroutines' equivalent von Mises stress-strain curves are nearly superimposed on top of one another for most of the flow stress curves. Furthermore, the flow stress curves and the distribution of equivalent strain and von Mises stress given by the models with 0.25 mm and 0.5 mm mesh sizes are indistinguishable. The values of equivalent plastic strain and von Mises stress at various sections of the model were also noted to correspond for these two mesh sizes. Given that more computation resources (space and time) are required as mesh density increases, the mesh of size of 0.5 mm was selected for simulation of the necking of a cylindrical numerical model at different strain rates and temperatures of 298 K and 773 K.
The results obtained from the last elements near the point y 2 (shown in Figure 16) were compared with the analytical solutions obtained from Equation (1) at different strain rates. The flow stress curves obtained from numerical simulation and the analytical constitutive model for sample types C, D, and E, in this case, are summarized in Figure 20.
Similar to the result from the single element test reported earlier on, the curves for the analytical solution and those from the VUMAT and VUHARD subroutines are nearly superimposed on top of one another at temperatures of 298 K and 773 K and respective strain rates, as seen in Figure 20. That the curves in this figure are nearly overlapping indicates that both subroutines are accurate in modelling a multiple element cylinder.
The foregoing discussion suggests that both subroutines implemented in this study can reliably be used to simulate the mechanical properties of various microstructures of additively manufactured Ti6Al4V parts. However, it has also been noted that the VUMAT subroutine is slower and requires more computational resources compared to the VUHARD subroutine. This is a clear indication that the different integration schemes employed by these subroutines influence the stable time increment.

Conclusions
The implementation and verification processes of the VUMAT and VUHARD subroutines for numerical simulations of plastic deformation of various forms of DMLS Ti6Al4V (ELI) were presented and discussed in this article. The following conclusions are deduced from the research:

Conclusions
The implementation and verification processes of the VUMAT and VUHARD subroutines for numerical simulations of plastic deformation of various forms of DMLS Ti6Al4V (ELI) were presented and discussed in this article. The following conclusions are deduced from the research:

•
Using single element numerical models and at a test velocity of 4 m/s, both VUMAT and VUHARD subroutines showed fairly similar values of average plastic strain rate of approximately 350 s −1 .

•
The equivalent plastic strain and von Mises stress generated from the two subroutines for the single element numerical model also coincide for a large part of the simulation time.

•
The single element numerical models developed in this study and simulated using the developed VUMAT and VUHARD subroutines were shown to be sensitive to the critical model parameters of initial dislocation density (ρ o ), test temperature (T), and strain hardening and the dynamic recovery parameters h and k 2 , respectively.

•
The results of simulation of samples C, D, and E using single element numerical models based on the VUMAT and VUHARD subroutines and those derived from the analytical solutions at the same conditions of strain rate and temperature were nearly overlapping over most of plastic strain range. The equivalent plastic strain and von Mises stress profiles for multiple element numerical models exhibited three typical zones which varied in range for different mesh sizes.

•
The values of equivalent von Mises stress obtained in VUHARD and VUMAT subroutines varied in the first and second zones for different mesh sizes. The small variations seen in these results were attributed to the different integration schemes employed by these subroutines.

•
The results obtained from both subroutines for the models with 0.25 mm and 0.5 mm mesh sizes were indistinguishable. • Comparisons between the results obtained for the multiple element numerical modelling of samples C, D, and E and the curves obtained from analytical solutions at the same conditions of temperature and strain rate were nearly superimposed.
It is recommended that future research aims at carrying out complex dynamic simulations of experimental tests, such as split-Hopkinson bar pressure tests, the Taylor impact test, and blast loading of plates using the developed VUMAT and VUHARD subroutines. The microstructure-sensitive model implemented as a material subroutine in this study should be extended to include some form of failure criterion or damage model.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/ma15010372/s1, S1: The developed code for VUHARD subroutine, S2: The developed code for VUMAT subroutine.  For the isotropic von Mises yield criterion, the yield surface can be expressed in terms of the second deviatoric stress invariant J 2 , and thus the yield function can be written as From the von Mises failure criterion (distortion energy theory) [34,35] and the second deviatoric stress invariant J 2 , Equation (A7) can be written as where σ y is the yield stress or the value of stress of the last increment. This scalar quantity of stress obtained from the deviatoric stress is commonly referred to as the von Mises equivalent stress, q, computed in most finite element analysis packages such as ABAQUS viz [27]: Likewise, the scalar quantity of strain from the deviatoric strain, commonly referred to as equivalent strain, is expressed as [27] ε eq = 2/3 e ij e ij (A10) When σ y − √ q > 0, the state of the material is no longer elastic, and therefore a plastic correction must be performed. In other words, the equivalent increment in inelastic strain must be calculated to return the isotropic stress state to a von Mises yield surface. The radial return mapping algorithm introduced by Wilkins [37] is the common method that uses the von Mises yield surface to define the increment of plastic strain necessary to return the stress state to the yield surface after an increment of trial stress (S t ij ).

Appendix A.2. The Radial Return Algorithms
It is a very rare situation for most materials to behave perfectly plastic, the state of which makes the yield surface expand during plastic deformation. In radial return mapping, the stress is updated with an assumption that the response is elastic, and if it is beyond the yield surface, the stress is projected onto the closest point on the yield surface. In elastoplastic deformation, the total change in strain (∆ε) can be expressed as the sum of the elastic (∆ε e ) and plastic strain (∆ε p ) changes: If the change in strain tensor (∆ε ij ) has been derived for the Gauss integration point from the current iteration, then the trial stress S t ij can be obtained as ij is the old trial stress obtained from previous iteration, and µ is the shear modulus of the material. The actual stress due to plastic deformation, obtained at the end of the increment, is not the trial stress but rather the one returned to the yield surface, as seen in Figure A1. From Figure A1, the actual deviatoric stress can be expressed as The correction, also referred to as return stress ( ), is expressed as The return stress ( ) leads to the name "radial return" method. The implicit integration scheme [34] is presented in Appendix B.
where ∂ p is the equivalent plastic strain increment, and q (2) is the new equivalent von Mises stress. Equations (A15) and (A17) can be combined into Thus, the new deviatoric stress is given by The new equivalent von Mises stress can be obtained from Equation (A19), thus making Equation (A8) of the appendix Since, at the yielding surface, q (2) = σ (2) y (where σ (2) y is the new stress computed from the constitutive equation), Equation (A21) is a non-linear yield function, and the increment of plastic strain ∂ p between states 1 (old state) and 2 (new state) must be computed. The iteration methods of finding solutions of non-linear functions can be used to do this. Alternatively, a noniteration method can be used to solve for ∂ p [27] if the plastic flow stress occurring between the old and new state can be considered. The flow stress at the new state can be expressed explicitly using the flow stress at the old state as σ (2) where H denotes the rate of change in flow stress with plastic strain. Equations (A21) and (A22) can be combined to solve for ∂ p , the new deviatoric stress, and the stress tensor as Once the increase in equivalent plastic strain is computed, the proportionality factor σ (1) y /σ t e is calculated using Equation (A24), and the trial deviatoric stresses are scaled back to the yield surface using Equation (A25). The new stress tensor is then updated by adding the scaled trial deviatoric stresses to the hydrostatic stress in Equation (A26).